function x_next_Nxzz = gen_law_of_motion_x(parms)

    NN = parms.NN;
    Nx = parms.Nx;
    Nz = parms.Nz;
    
    % dynamics for x' = x'(N,lambda,z,z')
    x_Nxz = permute(repmat(parms.grid_x(:),[1 NN Nz]),[2 1 3]);
    sigma_x_Nxz = permute(repmat(parms.sigma_x,[1 1 Nx]),[1 3 2]);
    x_next_Nxzz = zeros(NN,Nx,Nz,Nz);
    for iz_now = 1:Nz
        expected_znext = parms.StateTransitionProbs(iz_now,:)*parms.grid_z(:);
        if parms.OPTION_x_standardize_shock
            std_znext = sqrt(parms.StateTransitionProbs(iz_now,:)*(parms.grid_z(:).^2) - expected_znext^2);
        else
            std_znext = 1;
        end
        for iz_next = 1:Nz
            x_next_Nxzz(:,:,iz_now,iz_next) = (1-parms.rho_x)*parms.x_bar + parms.rho_x*x_Nxz(:,:,iz_now) ...
                + sigma_x_Nxz(:,:,iz_now)*(parms.grid_z(iz_next) - expected_znext)/std_znext;
        end
    end
    x_next_Nxzz = reshape(x_next_Nxzz,[NN*Nx*Nz,Nz]); % ((N,lambda,z),z') format
    x_next_Nxzz(x_next_Nxzz<parms.grid_x(1)) = parms.grid_x(1);
    x_next_Nxzz(x_next_Nxzz>parms.grid_x(end)) = parms.grid_x(end);

end